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Abstract 

Biological rhythms are generated by pacemaker organs, such as the heart pacemaker organ (the 
sinoatrial node) and the master clock of the circadian rhythms (the suprachiasmatic nucleus), 
which are composed of a network of autonomously oscillatory cells. Such biological rhythms 
have notable periodicity despite the internal and external noise present in each cell. Previous 
experimental studies indicate that the regularity of oscillatory dynamics is enhanced when noisy 
oscillators interact and become synchronized. This effect, called the collective enhancement of 
temporal precision, has been studied theoretically using particular assumptions. In this study, 
we propose a general theoretical framework that enables us to understand the dependence of 
temporal precision on network parameters including size, connectivity, and coupling intensity; 
this effect has been poorly understood to date. Our framework is based on a phase oscillator 
model that is applicable to general oscillator networks with any coupling mechanism if coupling 
and noise are sufficiently weak. In particular, we can manage general directed and weighted 
networks. We quantify the precision of the activity of a single cell and the mean activity of an 
arbitrary subset of cells. We find that, in general undirected networks, the standard deviation 
of cycle-to-cycle periods scales with the system size as 1/ \/A, but only up to a certain system 
size N* that depends on network parameters. Enhancement of temporal precision is ineffective 
when N > N* . We also reveal the advantage of long-range interactions among cells to temporal 
precision. 

Author Summary 

Various endogenous biological rhythms in our body such as heartbeats and sleep- waking cycles of 
about 24-hour period, the so-called circadian rhythm, function in our body. Unexpectedly, these 
rhythms maintain time regularly. For example, the daily onset of activity in mice has a standard 
deviation of a few minutes even in the absence of environmental information. These biological 
rhythms are generated by pacemaker organs composed of a network of autonomously oscillatory 
cells. How do biological cells generate highly precise rhythms despite internal and external 
noise present in each cell? We know, experimentally, that an isolated cell cannot generate such 
precise oscillation, but a network of coupled cells can. Regularity in oscillations increases with 
the number of cells that constitute the network. This effect is called the collective enhancement of 
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temporal precision. In this study, we present a new theory for quantifying temporal precision in 
terms of network parameters including the number of cells, connectivity, and coupling strength. 
Our main finding is that the collective enhancement is ineffective beyond a certain cell number, 
and this number increases with coupling strength among cells. Our theory provides a useful tool 
for inferring the properties of cell networks. 

Introduction 

Biological rhythms such as heartbeats and sleep- waking cycles are essential in living organisms. 
Many biological rhythms are generated by pacemaker organs composed of autonomously rhyth- 
mic cells. For example, the heart pacemaker (i.e., the sinoatrial node) is the source of electric 
waves propagating from within the heart, which cause the contraction of cardiac cells [T] The 
suprachiasmatic nucleus (SCN), which is a network of clock cells located in the brain, orches- 
trates the circadian (i.e., approximately 24 h) activity of the entire body. Each clock cell has 
a circadian rhythm in its electric activity owing to the gene regulatory network within the cell, 
and a population of clock cells synchronizes its activity through neural interactions [2]. The 
medullary pacemaker nucleus in electric fish is the pacemaker for the electric discharges emit- 
ted by electric fish, which are used for object detection and communication with other electric 
fish 0. 

Cell dynamics involve fluctuations resulting from various types of internal and external noise. 
However, oscillations in pacemaker organs such as the sinoatrial node in the heart, the SCN, and 
the medullary pacemaker nucleus in electric fish are highly precise. For example, the daily onset 
of activity in certain mammals and birds has a standard deviation (SD) of a few minutes even in 
the absence of environmental information [5]. In addition, the electric organ discharge pattern 
in certain electric fish has a standard deviation of as little as 0.02% of the average period (5]. 

Experiments by Clay and DeHaan provided an important clue for understanding the mech- 
anisms underlying precise oscillations was provided by [6]. They prepared clusters of cultivated 
cardiac cells, ranging in size from 1 to ~100, and observed the beatings of individual cells. 
They found that the SD of inter-beat intervals decreases with the number of component cells 
in the cluster (N) roughly as SD oc Therefore, precision in individual cell oscillations is 

enhanced as the number of cells increases. Note that this scaling, which is reminiscent of the 
central limit theorem, is not at all trivial. This is because oscillators are synchronized and thus 
strongly correlated, while the central limit theorem is applicable to an ensemble of independent 
elements. 

The decrease in SD as N increases, the so-called collective enhancement of temporal precision, 
has attracted considerable attention [4HT7]. There is a large body of experimental [5l[6l[8l[9], 
numerical |10H13j . and analytical |14H17j studies. Theoretically, it has been shown thatthe 
average activity of all oscillators on the all-to-all network (i.e., the complete graph) obeys SD oc 
l/\/iV [l31[T5]. However, most analytical studies are based on rather strong assumptions about 
coupling topology (e.g., all-to-all) or coupling mechanism (e.g., gap-junction type). Moreover, 
little is known about temporal precision in single cell activity or ensemble activity for a subset 
of cells in an entire network. Note that in the experiments by Clay and DeHaan, the behavior 
SD oc \/\fN was found for single cells and not for the entire network [B]. 
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In this paper, we propose a general theoretical framework that enables us to understand 
the dependence of temporal precision on network parameters, including size, connectivity, and 
coupling intensity. Our framework, based on a phase oscillator model, allows us to handle 
directed and weighted networks, various coupling mechanisms, and temporal precision in the 
activity of single cells and arbitrary subsets of cells. 

We begin by describing the numerical results for two biological pacemaker models: network 
of the FitzHugh-Nagumo oscillators and that of circadian oscillators. These models have distinct 
oscillation and coupling mechanisms. For different networks including all-to-all coupling, lattices 
with nearest-neighbor coupling, and the random graph, we observe that there is a common 
dependence of temporal precision on network size A^. The SD of cycle-to-cycle periods decreases 
as in small networks, but approaches an asymptotic value as N increases. That is, there 

is a crossover. Then, we develop a theory for obtaining an explicit expression for the SD of the 
cycle-to-cycle period. In particular, we find the condition for the behavior SD oc \/\fN and 
the dependence of the crossover point A^* on network parameters. We also demonstrate the 
advantage of long-range interactions among cells to temporal precision. Finally, we discuss the 
implications of our theory. 

Results 

Numerical results 

First, we present the numerical results for two mathematical models describing biological oscil- 
lations (see Methods for the details of the models). We used FHN oscillators with gap-junction 
coupling as a model of oscillatory cardiac or neural cells. We also employed a previously pro- 
posed model for the SON (i.e., a population of circadian clock cells) [18], which is referred to as 
the SCN model. 

Waveforms and oscillation periods are regularized when oscillators are coupled 

Figures [Ha,c,e) and (b,d,f) present the waveforms obtained from the FHN and SCN models of 
different network sizes, respectively. The average cycle-to-cycle periods are depicted by dotted 
lines in each panel to illuminate the variations in cycle-to-cycle periods. The properties of each 
constituent cell were kept constant, while the connectivity between the cells is different. Typical 
waveforms of uncoupled cells {N = 2) are shown in Figs. [Ha,b). When the cells are coupled 
sufficiently strongly, the system synchronizes stably (Figs. mc,d)). Figures [l][c,d) indicate that 
waveforms in the presence of coupling are regularized as compared to the waveforms of isolated 
cells [Figs. [Ha,b)]. In particular, the variation in the cycle-to-cycle period decreases. When 
100 oscillators are coupled [Figs. [Il[e,f)], the variation appears to be even smaller. When cells 
are coupled, individual cell oscillations are not only synchronized but also regularized, and the 
oscillation appears to be more regular for a larger system size. 

There is a limit to the enhancement of temporal precision 

To quantify the dependence of temporal precision on network parameters, we measured the 
coefficient of variation (CV), which is the SD of the cycle-to-cycle period divided by the mean 



4 



period. A cycle-to-cycle period is defined by an interval At between two successive passages 
of an observed variable (xj) across a specified threshold value xth (Fig. [2]). We set Xth = 0.4 
and 2.0 for the FHN and SCN models, respectively. We discard At that is much smaller than a 
typical oscillation period to exclude noise-driven rapid threshold crossing. The CV is defined as 



where r and SD are the mean and the SD of a series of At, respectively. 

Here we investigate the FHN model on networks of different types and different sizes. We 
assume that the system is composed of identical cells subjected to weak noise. Figure[3)I^a) shows 
the CV of individual cell oscillations in the FHN model on the all-to-all network of size A^. The 
results for different coupling strength values, k, are plotted using different symbols. We find 
that 

(i) CV is proportional to for small N values for each k 

(ii) CV approaches a constant value for large N values for each k; i.e., there is a crossover 

(iii) the crossover point N* increases with k. 

We observe similar behavior for the square lattice and the random graph, as shown in 
Figs. Eljb) and (c), respectively. 

Temporal precision increases with N, while the level of synchrony remains constant 

A natural question is whether the enhanced synchronization induces the collective enhancement 
of temporal precision. To examine this possibility, we measured the distance 6 between the 
actual state and the in-phase state (see Methods for the definition of 6) for the all-to-all 
network. As shown in Fig. [3]^d), the level of synchrony is independent of N for each k value. We 
also confirmed that, in the FHN model on a square lattice, 6 even increases with N although the 
CV decreases (results not shown). Thus, the enhancement of temporal precision by an increase 
in N is not attributed to the improvement in synchronization. 

CV for ensemble activity has a larger crossover point 

In nature, rhythmic output from a pacemaker organ is usually generated by an ensemble of 
multiple cells. For example, rhythmic electroactivity propagating within the heart is thought to 
originate from cells on the surface of the sinoatrial node. The SCN consists of various neural 
populations, and each population forms a particular pattern of efferent projections to other 
parts of the brain [19]. This anatomical fact suggests that the SCN's output is generated by a 
combination of a subset of neurons rather than by the uniform average of the entire organ. 

Therefore, we investigated the CV of the ensemble activity of a subset of cells on the all-to-all 
network. The ensemble activity is defined by the average waveform of M (1 < M < N) cells: 




M 



(2) 



i=l 
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where the measured ensemble is assumed to consist of oscillators xi, . . ., xm- The cycle-to-cycle 
period and the CV for the ensemble activity are defined similarly to the case of single cell activity 
(Fig. [2]). In Fig. [U we present the CV measured for the average waveform with different values 
of M in the FHN model on the all-to-all network. For M smaller than N, properties (i)-(iii) 
listed above are preserved. In addition, we find that 

(iv) the crossover point N* increases with ensemble size M 

(v) for M = N, the CV is proportional to for any N; i.e., there is no crossover. 

We also confirmed that the same properties hold true for the FHN model on two-dimensional 
lattices and for the SCN model on the all-to-all network and the two-dimensional lattice. 

Results are qualitatively the same under strong noise and heterogeneity 

So far, we have assumed an ideal case: identical oscillators and weak noise. To simulate more 
realistic situations, we now consider networks composed of heterogeneous cells subjected to 
relatively strong noise. As examples, we measure the CV for the FHN model on the square 
lattice and for the SCN model on the all-to-all network (Fig. [5]). In the FHN model, we made 
one of the parameter values heterogeneous in order to obtain the distribution of natural periods 
of cells as Tj w 133 it 3 (mean it SD). In the SCN model, the time scales of the cells were made 
heterogeneous such that Ti ~ 23.4 it 1.2. The latter situation is consistent with the experimental 
observation by Honma et al. [20]. In all cases, we apply sufficiently strong coupling to ensure 
that the oscillators are well synchronized. Under this condition, as seen in Fig. [5l all properties 
(i)-(v) hold true. 

Theory 

We found, numerically, that properties (i)-(v) hold true in various situations. In the following, 
we develop a theory for relating temporal precision to network parameters by assuming weak 
coupling and weak noise. Under this assumption, a large class of oscillator systems including the 
models considered above are reduced to the phase model (see Methods and References [2T|[22]) 
given by 

N 

4)i{t) = Ui + K^Aijf{(^j - (^i) + \^ii{t), (3) 
i=i 

where (pi and (1 < i < N) are the phase and intrinsic frequency of the ith. oscillator, re- 
spectively; A = {Aij) is the weighted adjacency matrix with its element Aij equal to the 
intensity of the coupling from the jth to ith oscillators; k is the overall coupling intensity; 
/(•) is a 27r-periodic function; ft(t) is independent white Gaussian noise with E[^j(t)] = and 
E[^j(i)^j (t')] = 5ij5{t — t'), where E represents the expectation; and D is the strength of the 
noise. The adjacency matrix A is allowed to be asymmetric, weighted, and to possess negative 
components. Extension of the following results in the case of i, j-dependent coupling function 
fij{-) and i-dependent noise strength Di is straightforward. For clarity of the presentation, we 
focus on Eq. We assume that all the oscillators are synchronized in frequency; i.e., all the 
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oscillators have the actual frequency $7 owing to the effect of coupling. Synchronization usually 
occurs when coupling is sufficiently strong compared to noise and heterogeneity in Wj. 

One oscillation cycle corresponds to an increase in the phase by 27r. More precisely, the A;th 
cycle-to-cycle period of the ith oscillator is defined by Atj'^^ = ~t^t where tf'^ is the first 
passage time for (j)i{t) to exceed 2k'K (Fig. [6|). Because we assumed that all the oscillators are 
synchronized to il, the expected value of At^*^^ (r) is independent of i and is given as 

r - E[Atf = ^, (4) 

where the statistical averages are taken over different k values. The temporal precision of the 
ith oscillator is characterized by 



SDi = std[At,] = Vvar[At,] = ^ E[(Atf ^ - rf] . (5) 
The CV for the ith oscillator is equal to 

CV. ^ ^. (6) 

T 

To obtain the dependence of SDj on network parameters, we employ an approximation given 

by 

27r 

— std[Ati] std[A(j)i\, (7) 
r 

where A(f)i = (f>i{t + t) — (f)i{t) — 27r (Fig. [7|). For an isolated oscillator obeying = uji + ^/D^i{t), 
one immediately finds that var[A(^j] = Dri, where Tj = In/uji. When oscillators are coupled and 
synchronized with frequency 17, we write 

var[A(^i] = fXiDr. (8) 

We refer to fj,i as the scaling factor of the ith oscillator (Fig. [7|) . 

To obtain an expression for /ij, we assume that noise is sufficiently weak and linearize Eq. ^ 
around the synchronized state. The synchronized solution (/'f(t) (1 < "i < N) is represented as 

^l{t) = nt + i^i, (9) 

where and are the constants derived by setting (pi = il. and D = in Eq. ([3]); i.e., 

N 

By introducing a small deviation 

Oiit) = Mt) - 4>m, (11) 

we obtain 

N 

Oiit) = K J^«;,,(0,- - 9,) + VD^iit), (12) 
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where Wij = Aijf'[ipj — ipi) is the effective couphng weight. For convenience, we rewrite Eq. (jl2p 
as 

N 

Oi{t) = -Kj2LijOj + ^ii{t): (13) 
i=i 

where L = {Lij) is the Jacobian matrix with its element given by 



Note that L has a zero eigenvalue with the corresponding right eigenvector u^^'^ = (1, . . . , 1)^ /^/N . 
Furthermore, because of the assumption of the stability of the synchronized state, the real parts 
of the other — 1 eigenvalues of L are positive, i.e., = Ai < ReA2 < • • • < RcAat. The as- 
sumption of the stability holds true when Wij > (1 < i,j < N) and the network described by 
the adjacency matrix {wij) is strongly connected |23fl25] . For more general cases, the stability 
condition is nontrivial. 

For in-phase synchrony (i.e., = for 1 < « < iV in Eq. ([9])), which occurs when the 
heterogeneity in the network and in individual oscillators is sufficiently small and/or the coupling 
is sufficiently strong, we obtain Wij cx Aij for 1 < i,j < N. In this case, L is the network 
Laplacian generalized for a directed and weighted network [26], given by 



L.^^i ^ r _ (15) 



-Aij for i / j 
T.i'^i for i = j 

Note that L is symmetric when the adjacency matrix A is symmetric. 

As shown in Methods, for any diagonalizable matrix L, we obtain /Uj = C^, where 

_ E[(g,(t + r)-g,(t))(g,-(t + r)-g,(0)] 

m,n (m+n>2) 

Here and v^^^ are, respectively, the right and left eigenvectors of L that satisfy 

the orthogonality and normalization conditions; i.e., Lu^^^ = \nU^"'\ v^'^^L = A„i)("\ and 

^(m)^{n) ^ j^^. y^^ ^ ^(m) . ^(n)^ 

For symmetric L, which is the case for in-phase synchrony on undirected networks, Eq. (jl6p 
becomes much simpler. Because all the eigenvalues are real, ii^") = v^^^~^ , Vmn = u^™-* • u^^^ = 
6mn for 1 < m,n < N, and Ei^i ^1"^ ''^^^^ ' ''^^"^ ~ fo^ n > 2, we obtain 

n=2 

Moreover, because of the normalization condition, Ei^i ''^i"^^!"^ ~ 1' mean of /ij over the 



entire network, (fi) = f^i/^^ is independent of the eigenvectors and is given by 

1 1 1 _ p-l^^nT 



iV iV ^ KA„r 

n=2 



8 



Crossover point N* increases with coupling strength k 

If the second term of Eq. (|18|) is negligible compared to the first term, we obtain (^u) ~ 1/-/V; i.e., 
the SD decreases proportionally to l/\fN. However, as increases, the second term becomes 
comparable at certain N* and even dominant for N ^ N* . If the eigenvalue spectrum converges 
to a certain density function q{\) as N ^ oo, we obtain 

{f^)^f^oo= / ^(A) ^ dX (Af^oo). (19) 

We later demonstrate the convergence for the all-to-all and ring networks. Spectra of finite 
dimensional lattices [27j, uncorrelated random graphs with arbitrary degree distributions |28j . 
and the small-world network with a fixed expected degree [29] also converge. By equating the 
first and second terms in Eq. ()18p . we estimate the crossover point as N* ~ l/^oo- Since /ioo 
monotonically decreases with increasing k, N* increases with n. 

Furthermore, if the second smallest eigenvalue A2 is nonvanishing in the limit — )• 00 (which 
is the case, for example, in the all-to-all network and various random networks including small- 
world networks |28ll29j ) and k is so large that e""'*'^'^ <^ 1, we obtain /ioo oc Then, the 
crossover point scales as 

N* oc K. (20) 



Crossover point N* is proportional to the size of a measured ensemble 

By assuming in-phase synchrony, we calculate the scaling factor of the noise reduction for the 
ensemble activity of an arbitrary set of oscillators. We rearrange the oscillator indices and write 
the ensemble activity as 

M 

X{t)=Y,C^X^{t), (21) 
i=l 

where Ci ^ is an arbitrary constant with the normalization condition '^ff^i Ci = 1. When the 
deviation 9i from in-phase synchrony (i.e., ipi = for 1 < i < iV in Eq. ([9])) is small for each 
oscillator, the phase of X(t) is approximated by 

M M 

$(t) = J^Ci'/'i(t) = f^t + I]Oei(t)- (22) 

i=l i=l 

Then, similar to the case of individual cell oscillations, we define the scaling factor for the 
ensemble activity as 

var[A$] = /i$L»r, (23) 



where A<I> = ^{t + r) - ^{t) - lir. We then obtain 
var|A5>| 



var[A^] _ ^ , , Emit + r) - g,(t))(g,-(t + r) - Ojit))] _ ^ , , ^ 

i,j=l i,j=l 

Henceforth, we assume Q = 1/M for 1 < i < M, as is the case in Figs. [3]and[5]^b). 
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There are notable properties for symmetric L (see Methods). When M = N (i.e., X{t) is 
the mean activity of the entire network) , we obtain 

^^^ = JJ^ (25) 

that is, there is no crossover. For M < N, fj,^ generally depends on the choice of M oscillators. 
However, if we randomly choose M oscillators out of N oscillators, where 1 ^ M ^ A^, we 
estimate 

1 1 A 1 - e-''^"^ 1 /Uoo 

n=2 



In this case, the lower bound of the SD is inversely proportional to v M and the crossover point 
increases as 

N* oc M. (27) 
As shown later, this estimation is asymptotically exact for the all-to-all network. 



Behavior can be violated even for small N values when L is asymmetric 

The behavior CV oc is obtained for N < N* when the Jacobian L is symmetric, which 

is the case when a network is undirected and the oscillators are synchronized in phase. We 
refer to this situation as "democratic" because symmetric L implies that the action and reaction 
between any two nodes are balanced. 

For asymmetric L, Eq. (jl6p implies that the SD at small N values decreases as \/Vii/N 
instead of 1/^/N. In [30], we analyzed the long-time diffusion property of Eq. ([3]) to obtain 
0-2 = limAt^oo var[6'i(i + At) - ei{t)]/{DAt) = Vu/N through a different technique. This 
previous result is consistent with that obtained in the present paper because o"^ corresponds to 
phase diffusion after infinitely many cycles, and the second term on the right-hand side of Eq. (jl6p 
vanishes with this limit. Furthermore, we showed in [30] that \/Vii/N is larger than or equal 
to l/\/iV for asymmetric L. For example, in directed scale-free networks, which is a strongly 
heterogeneous network, we obtained y^Vn/N oc A^"^ with < /3 < 1 /2; the effect of collective 
enhancement is significantly weaker. Moreover, the scaling ^/Vn/N = N~^^'^ can be violated 
even when a network is undirected. This is the case when the synchronized state is not in-phase 
but accompanies a wave pattern. Wave patterns arise when the network is spatially extended 
(such as Euclidian lattices) and the natural frequency is sufficiently heterogeneous [22l[31]. In 
this case, Vu/N decreases with N for small values but approaches a constant value for large 
values. Thus, strongly asymmetric connectivity and/or strong heterogeneity in the oscillator's 
properties can hamper the collective enhancement of temporal precision. 



Examples and numerical verification 

To demonstrate and numerically confirm our analytical results, we investigate the phase model 
(Eq. ([3D) on several networks. In numerical simulations, we set cj, = 1, f{(j)) = sIikJ), and 
y/D = 0.01 in Eq. ([3]). In the example networks, all the oscillators synchronize in phase in the 
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absence of noise. Thus, Wij = Aij and = a; for any coupling strength k and any N . Note that 
the dependence of the CV on k and is only through the SD because r = 2t: /u is constant. 
In the following, we show the values of the normalized CV, that is actual CV values divided 
by the CV of isolated oscillators, shown as ^/Dt/2tt. Our theory predicts that CVj « ^JJTi and 
CV^) ~ -y//x<j) . 

Two asymmetrically coupled elements. The first example is two asymmetrically coupled 
elements (A^ = 2): wi2 = p and W21 = 1 — p (Fig. [8]). In this case, we have Ai = 0, A2 = 
-l,uW = ^{l ,v^^^ = V2{p 1 - p),u^^'> = - p -p)"^, and 17(2) = -i=(l -1). By 

substituting them in Eqs. and for M = 2 and setting Ci = C2 = 1/2, we obtain 

m = ^2 = — +^1-— J^^, (28) 

where = 2p^ + 2(1 — p)^. For any k and r values, the best precision is obtained in the 
symmetric case {p = 0.5). Figure [8] suggests that the analytical and numerical results are in 
strong agreement. 

All-to-all coupling. The second example is all-to-all coupling; i.e., Wij = 1/N for 1 < i,j < N. 
The eigenvalues are given by = 1 (2 < n < A^). Because all the nodes are equivalent (i.e., 
permutation symmetry), we obtain /Xj = (fi). Then, from Eq. (jlSp . it follows that 

1 / 1 \ 1 - e~'^'^ 

We also obtain a concise form for /i$ (see Methods), given by 

1 /I 1 \ 1 - e~'^'^ 

We denote the CV value at = A^* by CV*. By equating the first and second terms on the 
right-hand side in Eq. (j3ip and assuming M <^ N and e~'^'^ ^ 1, we obtain 

A^* w ktM, CV* oc , ^ . (32) 

Figure [9] shows the analytical and numerical results. Note that in Figs. |4]and[Hb), the lower 
bounds are roughly proportional to as our theory predicts. 

Ring. The third example is the ring of size A'^, i.e., u)j_j+i = Wi^i-i = 1/2 for 1 < z < A^ and 
Wij = for j / i — 1, i + 1, as an example of spatially extended systems. For this network, we 
obtain 

A„ = l-™(?<!1^) (33) 

for 1 < ?i < A^. Because L is symmetric and the network has permutation symmetry, we obtain 
fii = {n) where {fi) is given by Eqs. (llSp and (jSSp . Figure [TOl shows the analytical and numerical 
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results. Although each cell is adjacent to just two cells for any > 3, there is a clear A^- 
dependence of the CV for individual cells. Temporal precision is not simply determined by local 
connectivity. 

The lower bound of the CV for the ring is considerably larger than that for the all-to-all 
network (Figs. El^a) and [TO]l . The reason for this is as follows. The Laplacian of the ring 
for a large A^ value has negligible eigenvalues (i.e., A„ for n ~ and n ~ A^ in Eq. (|33p ). 
and these eigenvalues significantly enlarge the second term of Eq. (Iisp . In contrast, there is a 
nonvanishing spectrum gap (i.e., the second smallest eigenvalue A2) in the all-to-all and various 
random networks |28tl29j. In the FHN model, we observed a similar difference between the 
cases of the square lattice (Fig. EJ^b)) and the all-to-all and random networks (Figs. [3]^a) and 
(c), respectively). This is also because the square lattice has negligible eigenvalues [27]. Such 
small eigenvalues are associated with slow synchronization of remote oscillators owing to a time 
lag in communication, and this property is shared by any spatially extended networks with 
local interaction. Therefore, spatial networks with only local interaction are disadvantageous to 
temporal precision. 

Small-world networks. By using a type of the Watts-Strogatz model |32p33j of fixed size A^, 
we demonstrate that a small fraction of long-range interactions added to the ring drastically 
improves temporal precision. We generate a network by adding pN bidirectional shortcuts 
sequentially to the ring, where p is the shortcut density. Under the condition that multiple 
links are avoided, the two endpoints of each shortcut are chosen from the A^ nodes with equal 
probability. The generated network is undirected. To maintain the total coupling strength 
independent of p, we set Wij = Wji = 1/(2 -|- 2p) for all links. The ring and all-to-all networks 
are obtained at p = and p = N/2 — 1, respectively. Figure [TT] shows the numerically obtained 
(CV) for each p, where (CV) = J^Zi CMi/N for a single realization of the network. The lines 
represent \/JJij obtained from Eq. ()18p . where we numerically computed the eigenvalues An for 
the generated network We set the coupling strength such that (/i) ^ 1/A^ (i-e., A^ ^ A^*) for 
the initial ring (p = 0). 

Figure [TT] indicates that temporal precision is considerably improved at p w 1, i.e., when 
0{N) shortcuts are added (the small- world regime). Moreover, the corresponding CV value is 
close to that of the all-to-all network, in which 0{N'^) "shortcuts" exist. As discussed above, 
there are small eigenvalues that hamper temporal precision in spatially extended networks. Such 
small eigenvalues do not exist in networks with a sufficient number of shortcuts because of rapid 
communication between any pair of oscillators. 

Mechanism of the crossover 

We demonstrated using various models that the crossover generally occurs in the collective 
enhancement. On the basis of our theory, the crossover can be interpreted as follows. When 
Jacobian L is symmetric, the SD for the mean phase decreases as l/\/iV for any A^ (Eq. (j25p ). 
When coupling strength k is infinite, the oscillators are completely synchronized in phase. Then, 
the phase of each oscillator is identical with the mean phase, and so is the SD, i.e., SDj oc l/\fN 
for any N. This behavior is expressed in the first term on the right-hand side of Eq. (jl7p . 
However, for finite k, individual oscillators' phases fluctuate around the mean phase because of 
the independent noise applied to the oscillators. Owing to this additional fluctuation, the SD for 
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individual oscillators is larger than that for the mean phase, as expressed in the second term on 
the right-hand side of Eq. (jl7p . Although the fluctuation in the mean phase vanishes with the 
limit — )• oo, it remains finite in individual oscillators. This is the origin of the lower bound. 

Discussion 

We found that the collective enhancement is ineffective for system size N above the crossover 
point A^*. We further showed that A^* increases with coupling strength (Eq. (j20p ). Therefore, as 
oscillators are more strongly coupled, the behavior CV oc 1/VN persists up to a larger A^ value. 
This is the case for different oscillation and coupling mechanisms, as demonstrated in the two 
biological models (the FHN and SCN models) and the phase oscillator model. Moreover, this 
behavior also holds true for different network connectivities, as demonstrated using the ring, the 
square lattice, the all-to-all network, and the random graph. 

Our theory is useful for inferring the magnitude of fluctuations in individual cells and the 
coupling strength between cells. Suppose that temporal precision in a pacemaker tissue that is 
genetically modified or subjected to a treatment (e.g., drug) is lower than that in an intact tissue. 
If the cells in the tissue are well synchronized in both cases, one may consider that the treatment 
affects the oscillation mechanism of individual cells. Our theory suggests another possibility: a 
decrease in the coupling strength, not the alteration in the oscillation property of individual cells, 
may be the reason for the reduced temporal precision (Figs. [3l|^a,b,c)). By observing reduced 
temporal precision only, we cannot distinguish these two possibilities. However, our theory makes 
it possible to individually quantify the effects of the treatment on the two properties if we can 
observe cell networks of different sizes. By observing temporal precision in small (i.e., A^ < A^*) 
tissues of different sizes, we can infer the magnitude of fluctuations in individual cells by fitting 
the law CV oc 1/^/N. Furthermore, by observing relatively large tissues and determining A^* 
values for different treatments (e.g., different days of cultivation, different concentrations of a 
drug, treated versus untreated), we can infer changes in the coupling strength induced by the 
treatment because A^* increases with the coupling strength (Eq. ()20p ). 

Our study also indicates that long-range interactions among cells are advantageous to tempo- 
ral precision. As demonstrated in Fig. [TT| the addition of shortcut links considerably decreases 
the CV. A similar result was reported in a previous numerical study using a more realistic model 
for the SCN jl3j . This result might underlie an evolutionary origin of dense fibers across the 

SCN dg. 

Our theoretical results provide an interpretation of previous experiments on cardiac and 
circadian oscillations. Kojima et al observed a decrease in the CV with increasing cell number 
in cultivated cardiac cells coupled via micro channels [9]. They showed that the CV decreases 
considerably with A^ for small A^ values (A^ = 1,2,3), while it is almost constant for A^ > 4. 
In contrast, in cultivated cardiac cells that are directly and tightly coupled to each other. Clay 
and DeHaan found that the reduction in the CV roughly obeys CV oc l/^fN up to A^ ~ 100 [6]. 
Although the cells are kept synchronized in both cases, the behavior of temporal precision is 
different. This discrepancy may be due to a difference in coupling strength. While the coupling 
was strong enough to guarantee synchrony in both cases, coupling in the latter experiments may 
be stronger than that in the former experiments, resulting in A^* ~ 4 and A^* > 100, respectively. 
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It would be of great interest to investigate systematically how the crossover point increases with 
coupling strength, possibly controlled by the width of the micro channel implemented in the 
former experiments [9]. 

Collective enhancement has been examined experimentally in circadian oscillation as well. 
Herzog et al measured temporal precision in SCN cells [8]. There, individual cell oscillations in 
both synchronized and unsynchronized cases were observed in slice cultures of SCN and dispersed 
SCN cells, respectively. They found that the SD in the former (0.42 h) was approximately 
five times smaller than that in the latter (2.1 h), and argued that the collective enhancement 
of temporal precision occurs in synchronized cells. They further speculated that, under the 
assumption SD ~ l/\/iV, only 25 cells out of the order of 10^ cells composing the SCN are 
involved in the collective enhancement of temporal precision in the explant SCN. 

We interpret this experimental result as follows. In the SCN, a wave pattern is observed 
[3^135] . As indicated above as well as in our previous paper [30] , the law SD oc 1/ \/iV is violated 
in the presence of a wave pattern even if the coupling is sufficiently strong. Roughly speaking, 
the reason for this is that only the cells forming the source of the wave pattern can contribute 
to the collective enhancement of temporal precision, and other cells simply obey those cells |30j . 
The number of cells forming the source might be of the order of 25. Cells located downstream 
of the wave may contribute to functions other than temporal precision. 

Our theory is widely applicable to frequency-synchronized oscillators with weak noise and 
weak coupling. Our theory can also apply to the case of the coexistence of multiple coupling 
mechanisms, only by replacing coupling function / by fij in the phase model (Eq. ([3])). Although 
the phase model is not justified when the assumption of weak noise and weak coupling is violated, 
we have numerically confirmed that our main finding, i.e., the properties (i)-(v), are preserved 
in the case of strong coupling and strong noise (Fig. [5]) . We thus expect that our theory, based 
on the phase model, captures the essence of the collective enhancement of temporal precision. 

Methods 

Model equations for biological pacemaker systems 

We consider two systems — the FHN model and the SCN model representing the cardiac pace- 
maker organ and the circadian master clock, respectively. 

The FHN model has been extensively used as a model of neurons and cardiac cells [36]. Our 
FHN model is given by 

TV 

Xi){xi - I) -yi + pCi{t) + K'^Aij{xj - Xi), (34a) 

i=i 

byi + c), (34b) 

where a,b,c,e are the model parameters, p is the noise strength, £,i{t) is white Gaussian noise 
with E[^j(t)] = and E[^j(i)^j(t')] = 5ij6{t — t'). We chose parameter values such that each unit 
is autonomous oscillator. In Figs. [H [3]and[4l we set a = 0.1, e = 0.01,6 = 0.5, and c = 0.05. 
In Fig. [5]^a), we replace c with q = 0.1 + 0.02z/j (1 < i < N), where z/j is a random variable 



— — = xAa 
dt ^ 



dm 

dt 



e[Xi 
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independently taken from the Gaussian distribution with zero mean and unit variance. We 
varied the noise strength and coupling strength, as specified in the figures and their captions. 
The distance 5 from the in-phase state is defined as 



1 ^ 



(35) 



where x = X^^i Xi/N . 



As the SON model, we employed a previously proposed model [TH], given by 

at V + Zf K2+Xi Kc + KFiJ 

dyi_rr.(u.^. TA Vi (36b) 



dt \ K4 + y. 



dt \ Kq + z- 



dr. 









hVi - 




k-jXi - 


TV 




E 









+ P&\ (36c) 



dt \ Ka + ri 



(36e) 



where Vi = 6.8355, n = 5.6645, i^i = 2.7266,^2 = 0.2910, A:3 = 0.1177, V4 = 1.0841,1^4 = 
8.1343, fcs = 0.3352,1/6 = 4.6645,1^6 = 9.9849, /cr = 0.2282,1/8 = 3.5216, i^s = 7.4519,14 = 
6.7924, Kc = 4.8283, k = 12.0, and V2 = 12.0. All parameter values except for k and V2 are 
taken from [18]. Time constant Tj is introduced to express heterogeneity in the oscillation period. 
We set Tj = 1 in Fig.[TJ In Fig.[5)^b), Tj = l + 0.05i^j with Vi independently obeying the Gaussian 

distribution with zero mean and unit variance. The functions s!f\t) {C, = x,y,z,r) represent 
white Gaussian noise processes with E[^P(t)] = and {t)^^'^\t')] = 5ij5c_r^5{t - t'). The 

noise strength p and coupling strength k are specified in the figures and their captions. 

In both models, we applied sufficiently strong coupling to ensure that the oscillators were 
synchronized nearly in phase. When we computed the CV, we assumed random initial conditions 
and measured a sufficiently large number of cycle-to-cycle periods after the transient. 



Networks 

The all-to-all network used in Figs. [l]^b,d,f), [3l^a,d), HI [Sjb), and [9] is defined by Aij = 1/N for 
^ ^ h 3 ^ N . The one-dimensional lattice with an open boundary condition used in Fig. HJe) 
is defined by Aij = 1/2 for 1 < i < and l<j = iibl<A^, and Aij = otherwise. The 
ring used in Fig. [10] is the same as the one-dimensional lattice except that we impose a periodic 
boundary condition Ai^n = ^TV,i = 1/2- The square lattice with an open boundary condition 
used in Figs.[3l^b) and[5l^a) is defined by Aij = 1/4 with cell j adjacent to i for 1 < i,j < ^/N and 
Aij = otherwise. The undirected random graph used in Fig. Ws) is the Erdos-Renyi random 
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graph, where Aij = Aji = 1/8 for I < i < j < N with probabihty p = 8/N and Aij = Aji = 
otherwise. We set hnk weights such that the summed weight of the hnks per node is independent 
of N; i.e., for 1 < i < N , X^jLi ^ij = 1 in the ring and all-to-all network and X^jLi ^ij ~ 1 in 
the other networks including the Watts-Strogataz model used in Fig. [TTl 



Phase description 

A large class of oscillator systems including the FHN and SCN models (Eqs. (|34p and (j36|) ) are 
reduced to phase models if the coupling and noise are sufficiently weak |2HI22|. The concept 
behind the reduction is as follows. We denote an element of the state variable of the ith 
oscillator by Xi{t). When unperturbed, the oscillator portrays a one-dimensional closed orbit 
after transient so that Xi{t) = Xi{t + 2-k / uJi) , where cjj is the intrinsic frequency. We define the 
phase (pi by Xi{t) = Xi{(j)i/uji); that is, the phase increases linearly with time in the unperturbed 
oscillator. For convenience, we denote the unperturbed orbit by xi'Pi) = ^ii^Pi/^i)- Although 
the trajectory deviates from the closed orbit when the oscillator is weakly perturbed, it is still 
possible to parameterize a trajectory of an oscillator by only the phase and describe the dynamics 
of coupled oscillators in terms of the phases only |2HI22j. The resulting equation is given by 
Eq. ([3]). Because of the assumption of weak perturbation, Xi{t) is approximated by that of the 
unperturbed orbit, i.e., 

x^{t)^x{Mt))■ (37) 
Therefore, the first passage time problem for Xi{t) is approximated by that for (j)i{t). 

Calculation of Eq. ^ 

Our linearized equation is given by Eq. (|13p. which is reproduced as 

N 

9i{t) = -^Yl + ^ei(i), (38) 

where Oi {1 < i < N) is the deviation from the synchronized state, k > is the coupling strength, 
L = (Lij) is a diagonalizable matrix, and (,i{t) is white Gaussian noise with 

E[ei(t)]=0, E[m^j{s)]=6ij5{t-s). (39) 

From the assumption of the stability of frequency synchronization, we have 

= Ai < ReA2 < ReAs < ••• < Aat. (40) 

The right and left eigenvectors of L corresponding to An are denoted by u^^^ = (u^"^) and 
= {v^^'^), respectively; i.e., 

Ln(") = A„w("), (41a) 

= XnV^''^ (41b) 
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with the normahzation and orthogonahty conditions 
Using these eigenvectors, we decompose 6i{t) as 



TV 



where (pm{t) is given by 



m=l 



N 



i=l 



By taking the time derivative of Eq. ()44p and using Eqs. (j38p . (j42p and (j43p . we obtain 
where 



N 



1=1 



Equation (p9]) yields (r/m(i)) = 0. We also have 



E[r/^(t)?7n(s)] = E 



i=l 



TV 



TV 



DVmnS{t - s), 



6{t - s) 



where 



N 



(m) (n) 



i=l 



Now we derive Cij given in Eq. ()16p . The definition of Cjj is 



(42) 
(43) 

(44) 

(45) 
(46) 



(47) 



(48) 



(49) 
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By substituting Eq. (jl3]) in Eq. (jl9j) . we obtain 

N 

Dt 



1 ^ 

= ^Y. ^^"^^f [(v9„,(t + r) - <^^;„(t))((^„(t + r) - . (50) 

m,n=l 

The solution to Eq. (|45p is formally written as 

^^{t) = e-^^'-VmW + f ds e-^^-^'^'-'^M- (51) 

Jo 

Using Eq. (|5T]l . we obtain 

(^i(t + T) -(^i(t) = / dsi]i{s), (52) 



dsi dsaE [771(51)771(52)] 

rt+T rt+T 

= DVii / dsi / (is2 (5(si - S2) 
Jt Jt 

= DVuT. (53) 
To evaluate the terms on the right-hand side of Eq. (j50p for tti + 7i > 2, we first calculate 

Bmn{tl,t2) = f'dsi r ds2e-^^^-^'^-'^'^-^^-^''-'''>E[r,MVn{s2)] 

g-KAm(tl-s)-KA„(t2-s) 

(54) 





rt2 


/ dsi 


/ ds2 e 


Jo J 


'0 






/ dsi 


/ ds2 e 


Jo J 


'0 




min(ti,t2) 


DVmn [ 




Jo 






-KAm{tl — 


DVrr.„ ^ 





K{Xm + An) 

s=min(ti ,£2) 

We consider the limit t — )■ 00 in Eq. ()54p because we are concerned with a stationary process. 
Using Eq. ()54p . we obtain 

E [(99„,(t + r) - ¥P„^(t))((^„(t + r) - (^^(t))] 

= Bmn{t + r, i + r) - -B„„(t + r, t) - Bmn{t, t + T)+ Bmn{t, t) 

7T (t^Oo). (55) 



By combining Eqs. ()50p . (j53p and (|55p . we obtain 



C.^-FiiT., 7.^. + 2^ ^ VmnU, . (56) 

m,n (m+n>2) 

In Eq. ([H]), we set lif ^ = l/\/]V for 1 < i < iV. 
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Scaling factor for the ensemble activity 

In this section, we derive used in Eqs. ([25]) and ([25]) . By substituting Eq. ([57]) in Eq. ([^ . 
we express the ensemble activity X{t) in terms of the phases as 

M M 

X(t) = j;C.^,(t)«5^CiX(</'i(t))- (57) 

1=1 i=l 

For in-phase synchrony (i.e., ipi = 0) and smah deviation 9i, we can further approximate X(t) 
to 

M 

Xit) « xm + E Ca'(f^i)^^(i) « Xim), (58) 
where x'(</') = d,x{(t>) / d(p and <I> is the mean phase of the ensemble, given by 

M 

^it) = nt + Y,CMt)- (59) 

1=1 

Thus, similar to the case of individual cell oscillations, the cycle-to-cycle period for the ensemble 
activity X{t) is approximated by the cycle-to-cycle period At^^ for the mean phase $(t). We 
further employ the following approximation (Fig. [7j) 

2tt 

— std[At$] ^ std[A$], (60) 
r 

where A$ = ^>(t + r) — $(i) — 2tt. We define the scaling factor for the ensemble activity as 

var[A$] = fi^Dr. (61) 

We then obtain 

var[Acl>] ^ ^ ^ Emit + r) - 0.(t))(g,-(t + r) - g,-(t))] .... 

i,j=l i,j=l 



where Qj is given by Eq. (|T6|) . 

We consider the case of symmetric L and = 1/M for 1 < i < Af. Substituting Eq. (I17p in 
Eq. (j62p . we obtain 

71=2 1,^ = 1 

For M = N, nj"^ = ti^^) • = for 2 < n < TV (orthogonality) leads to 

= ^> (64) 
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that is, there is no crossover. For M < N, Eq. (j63p implies that /i$ depends on the choice 
of M oscillators. When we randomly choose M out of oscillators, where 1 <C M <C A'^, 
the dependence of /i$ on M is estimated as follows. The orthogonality and normalization, 
respectively, imply 

N N 
J- (n) „ J- (n) (n) ^ 

Therefore, the distribution of u["'^ (1 < ^ < N) has the mean of and variance of 1/A^. We 
randomly choose M N) elements and assume that they are independent random numbers 
with the same mean and variance. Then, we apply the central limit theorem for M S> 1 to 
obtain 

M M M 
i=l j = l 1=1 

By substituting Eq. (f66]l in the right-hand side of Eq. (1631) . we obtain Eq. (f26]l . 



Calculation of Eq. §3) 

It is convenient to choose the eigenvectors = {u'f'^) for 2 < ?i < as u[^^ = Xj^JvP' — n for 

1 < f < n — 1, ni"^ = (1 — n)/ \J v? — n and = for n < i < N . Then, the following property 
holds: 

I ^ ^ ^ ( for 2 < n < M, 

M^^i ^ ] ^ for M + 1 < n < iV. ^^^"^ 

Substitution of Eq. (f67l) and the eigenvalues A„ = 1 (2 < n < A^) in Eq. (|63|) results in 

N / , M 



n=2 ^ 

_ 1 1 - e ''^ 



n=2 \"' i=l 




N KT ^ vP- — n 

n=Af+l 

j_ 1 - e-"^^ A / 1 \_ 
N KT ^-^ \n — 1 n 

1 /I 1 \ 1 - e"''^ 



Af VM Af/ 



(68) 
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Figure Legends 



FHN model 



SCN model 






\ f 




1 f 








1/ 




1^ 




o o 



oo 



ooo 
ooo 




Figure 1. Waveforms obtained from biological oscillator models. We present the time 
series of Xi{t) of (a,b) two isolated cells (k = 0, N = 2), (c,d) two coupled cells («; > 0, = 2), 
and (e,f) two cells in 100 coupled cells {n > 0, N = 100) in (a,c,e) the FHN model and (b,d,f) 
the SCN model. In (e), we employ the one-dimensional lattice with an open boundary 
condition and show the waveforms of two neighboring cells. In (f), we employ the all-to-all 
network (Aij = 1/N for 1 < i,j < N). We set p = 0.1 in all panels and (a,b) k = 0, (c,e) 
K = 2, and (d,f) k = 1. 
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Figure 2. Schematic illustration of the concept of cycle-to-cycle period. 
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Figure 3. CV for single cell oscillations and synchronization distance 6 in the FHN 
model. (a,b,c) CV values for single cell oscillations on (a) the all-to-all network, (b) the 
square lattice, and (c) the undirected random graph of size A'^. (d) Distance 6 from in-phase 
synchrony in the all-to-all network. In (a) and (d), we set Aij = 1/N for 1 < i,j < N. In (b), 
the CV of the oscillator at the center of the square lattice with an open boundary condition is 
presented. We set Aij = 1/4 with cell j adjacent to cell i and Aij = otherwise. In (c), the 
CV value at given k and N values is defined as (CVi) = CVi/N for a single realization of 

the network. We set Aij = Aji = 1/8 with probability p = 8/N {I < i < j < N), and Aij = 
otherwise. The lines are guides to the eyes. We considered identical cells and weak noise 
(p = 0.01). The average period r is almost constant (r ~ 177) irrespective of N and k. 
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Figure 4. CV for ensemble activity of M cells in the FHN model on the all-to-all 
network. Parameter values are the same as in Fig. [3l except p = 0.0256 and K = 0.2. The 
hne is a guide to the eyes. 
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Figure 5. CV for biological models composed of heterogeneous cells subjected to 
relatively strong noise, (a) CV for single cell oscillations in the FHN model on the square 
lattice. We set p = 0.09. (b) CV for the ensemble activity of M cells in the SCN model on the 
all-to-all network. We set p = 0.04 and K = 12. The all-to-all network and the square lattice 
are the same as those in Figs. [3ja) and (b), respectively. The lines are guides to the eyes. 



29 



(a) uncoupled 

std[A^,] 



2n 




(b) coupled 




std[A(/),] 



time 




std[A(/),] 



time 



2n 



std[A/J = std[A(().] = ylDi 



2tz I 

— std[A?.] = std[A(l).] = y)J../>T 



Figure 7. Schematic illustration of our approximation to the cycle-to-cycle 
variation. Green trajectories represent different realizations of the phase (pi{t) of a single 
oscillator in (a) uncoupled and (b) coupled cases, where we set 0i(O) = 0. The red curves on 
the top of each panel represent distribution function Pt{t) obtained from the first passage time 
of 4'i{t) = Itx in different realizations. Our concern is its standard deviation, SDj = std[Atj]. 
The blue curves on the right of each panel represent distribution function P^{(j)i) obtained 
from different realizations of <j)i{T), where r is the mean period, and its standard deviation is 
denoted by std[A(^j]. We approximate std[Aij] using std[A0j]. Suppose that 
Pt{t)dt = Pff){4>i)d(j)i. On average, the phase crosses 27r with slope 2it/t, i.e., d(j)i/dt = Itt/t. 
We thus obtain Eq. d?]). For uncoupled oscillators (k = 0), our model corresponds to the 
Wiener process with a constant drift. In this case, Eq. ([7]) is exact, and we obtain 
(27r/r)std[Atj] = std[A(/)] = \/ Dt [37j. We also know that Eq. ([7]) is asymptotically exact in 
the one-dimensional Ornstein-Uhlenbeck process for weak noise [37]. For coupled oscillators 
(k > 0), however, our model (l3|) is a multivariate Ornstein-Uhlenbeck process when linearized. 
Even in this case, as is numerically confirmed in the examples shown in Figs. ISl llll Eq. ([7]) 
provides a suitable approximation. 
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Figure 8. Normalized CV versus coupling strength in asymmetrically coupled 
phase oscillators {N = 2). Presented is the normalized CV, i.e., CVj/(\/Dr/27r) {i = 1,2) 
and CV$/(\/Dr/27r) for two coupled phase oscillators with (a) p = 0.5 and (b) p = 0.2. 
Numerical results are shown by symbols. The solid and dotted lines represent the analytical 
results given by Eqs. ^ and ([29]), respectively Note that Vn = 2p^ + 2(1 - pf = I for 
p = 0.5 and Vu = 1.36 for p = 0.2. 
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Figure 9. Normalized CV in phase oscillators on the all-to-all network, (a) CV for 

individual cells (M = 1) for various k and N values, (b) CV for ensemble activity for various 
M values with k = 0.5. Symbols represent numerical data. Solid lines represent (a) y/JIi given 
by Eq. (130D and (b) ^//xi given by Eq. 
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Figure 10. Normalized CV for single cells in phase oscillators on the ring. Symbols 
represent numerical results. Solid lines represent given by Eq. (jlSp with Eq. (|33p . 
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Figure 11. Normalized CV for individual cells in phase oscillators on the variant 
of the Watts-Strogatz model. We present {C\i) /{V^ /2tt), where (CVi) = Eili CVi/A^. 
Lines represent y/JIi given by Eq. (fTTl) . We set N = 400. 



